function [ output_args ] = PRE5_GenerateAverageSubjects( input_args )
%PRE5_GENERATEAVERAGESUBJECTS Summary of this function goes here
%   Will generate the mean for all timefrequency files and create new
%   subjects to store this information in. Then it averages the results
%   into more new subjects.
%   The final result is a set of subjects (avg...) that each have an amb.
%   stim., amb.resp. and unamb.resp. datafile.
%
%   Note: This code contains information that is hardcoded, such as the
%   reference to the source subjects and the naming of the new subjects.

% Access the preloaded database..
global Database;

warning('PRE5_GenerateAverageSubjects: HARDCODED STUFF!!');

%% Create CWSphere
CreateNewMeanedSingleTrialSubject('CW','Sphere Ambiguous Response','CWSphere','Ambiguous Response');
CreateNewMeanedSingleTrialSubject('CW','Sphere Unambiguous Response','CWSphere','Unambiguous Response');
CreateNewMeanedSingleTrialSubject('CW','Sphere Unambiguous Stimulus','CWSphere','Unambiguous Stimulus');

%% Create CWFaceHouse
CreateNewMeanedSingleTrialSubject('CW','FaceHouse Ambiguous Response','CWFaceHouse','Ambiguous Response');
CreateNewMeanedSingleTrialSubject('CW','FaceHouse Unambiguous Response','CWFaceHouse','Unambiguous Response');
CreateNewMeanedSingleTrialSubject('CW','FaceHouse Unambiguous Stimulus','CWFaceHouse','Unambiguous Stimulus');

%% Create DSFaceHouse
CreateNewMeanedSingleTrialSubject('DS','FaceHouse Ambiguous Response','DSFaceHouse','Ambiguous Response');
CreateNewMeanedSingleTrialSubject('DS','FaceHouse Unambiguous Response','DSFaceHouse','Unambiguous Response');
CreateNewMeanedSingleTrialSubject('DS','FaceHouse Unambiguous Stimulus','DSFaceHouse','Unambiguous Stimulus');

%% Create JFSphere
CreateNewMeanedSingleTrialSubject('JF','Sphere Ambiguous Response','JFSphere','Ambiguous Response');
CreateNewMeanedSingleTrialSubject('JF','Sphere Unambiguous Response','JFSphere','Unambiguous Response');
CreateNewMeanedSingleTrialSubject('JF','Sphere Unambiguous Stimulus','JFSphere','Unambiguous Stimulus');

%% Create avg subjects
MergeMeanedSingleTrialSubjects({'CWSphere', 'CWFaceHouse'}, 'avgCW');
MergeMeanedSingleTrialSubjects({'CWSphere', 'JFSphere'}, 'avgSphere');
MergeMeanedSingleTrialSubjects({'CWFaceHouse', 'DSFaceHouse'}, 'avgFaceHouse');
MergeMeanedSingleTrialSubjects({'avgCW', 'JFSphere', 'DSFaceHouse'}, 'avg');

end


%% Support functions

function CreateNewMeanedSingleTrialSubject(sourceSubject, sourceCondition, destSubject, destCondition)
% Create a new subject based on a mean over epochs from the original
% subject. E.g.:
%    "CW, Sphere Ambiguous Response"  (39 epochs) --> "CWSphere, Ambiguous Response" (1 epochs)
%    "CW, FaceHouse Ambiguous Response"  (39 epochs) --> "CWFaceHouse, Ambiguous Response" (1 epochs)
% 
% Do this for all channels available for this subject/condition source.
%
% Goal:
% The new subjects (CWSphere, CWFaceHouse, JFSphere, DSFaceHouse) each have
% a Amb.Stim, Amb.Resp and Unamb.Resp. file for every channel.
%
% Note: the function writes its result in the timefrequency data folder and
% updates the timefrequency datatable accordingly.

% Setup
global Database
cFilename = strcmp(Database.timefrequency_header, 'Filename');
cSubject = strcmp(Database.timefrequency_header, 'Subject');
cChannel = strcmp(Database.timefrequency_header, 'Channel');
cChannelType = strcmp(Database.timefrequency_header, 'ChannelType');
cCondition = strcmp(Database.timefrequency_header, 'Condition');
cFrequencies = strcmp(Database.timefrequency_header, 'Frequencies');
cTimes = strcmp(Database.timefrequency_header, 'Times');
dtTimefrequency = Database.timefrequency;
% Select the correct rows
tfSubject = strcmp(dtTimefrequency(:,cSubject), sourceSubject);
tfCondition = strcmp(dtTimefrequency(:,cCondition), sourceCondition);
channels = unique(dtTimefrequency(tfSubject, cChannel));
% Loop over channels
for i = 1:length(channels)
    channel = channels{i};
    CurrentRow = tfSubject & tfCondition & strcmp(dtTimefrequency(:,cChannel), channel);
    % Load and transform the data
    filenameSource = dtTimefrequency{CurrentRow, cFilename};
    SOURCE = load(filenameSource);
    power = mean(SOURCE.power,3); % AVERAGE OVER TRIALS
    % Write/copy other properties
    nrOfEpochs = 1;
    subject = destSubject;
    channelType =  dtTimefrequency{CurrentRow, cChannelType};
    condition = destCondition;
    frequencies =  dtTimefrequency{CurrentRow, cFrequencies};
    times =  dtTimefrequency{CurrentRow, cTimes};
    % Write to file and update in database
    filenameDestination = [Database.Path.TimeFrequencyData subject ' ' channel ' ' condition ' ' num2str(min(frequencies)) '-' num2str(max(frequencies)) 'x' num2str(length(times)) '.mat']; %#ok<NOPRT>
    save(filenameDestination, 'power', 'frequencies', 'times',  'channel', 'channelType', 'subject', 'condition', 'nrOfEpochs', '-v7.3');
    UpdateOrAddTimeFrequencyFileInDatabase(filenameDestination, subject, channel, channelType, condition, frequencies, times, nrOfEpochs)
    % Celebrate the succesfulness
    disp(['- Created ' filenameDestination]);
end

end



function MergeMeanedSingleTrialSubjects(subjects, newSubject)
% Merge these subjects into a new subject.
% E.g.:
%   Merge({'CWSphere', 'CWFaceHouse'}, 'avgCW');
%
% The function will do this for all conditions and all channels found in
% both source subjects.
%
% Note: the function writes its result in the timefrequency data folder and
% updates the timefrequency datatable accordingly.

% Setup
global Database
cFilename = strcmp(Database.timefrequency_header, 'Filename');
cSubject = strcmp(Database.timefrequency_header, 'Subject');
cChannel = strcmp(Database.timefrequency_header, 'Channel');
cChannelType = strcmp(Database.timefrequency_header, 'ChannelType');
cCondition = strcmp(Database.timefrequency_header, 'Condition');
cFrequencies = strcmp(Database.timefrequency_header, 'Frequencies');
cTimes = strcmp(Database.timefrequency_header, 'Times');
dtTimefrequency = Database.timefrequency;

% Create the SELECT variables for subject
for sIdx = 1:length(subjects);
    tfsSubject{sIdx} = strcmp(dtTimefrequency(:,cSubject), subjects{sIdx});
    if ~any(tfsSubject{sIdx}); warning(['No data for subject ' subjects{sIdx} '.']); end;
end

% Find the conditions and channels that are common across the source subjects..
allChannels = [];
allConditions = [];
for sIdx = 1:length(subjects);
    allChannels{sIdx} = unique(dtTimefrequency(tfsSubject{sIdx}, cChannel));
    allConditions{sIdx} = unique(dtTimefrequency(tfsSubject{sIdx}, cCondition));
end
baseChannels = allChannels{1};
tfAllChannels = true(1, length(baseChannels));
for sacIdx = 1:length(baseChannels)
    for sIdx = 2:length(subjects)
        tfAllChannels(sacIdx) = tfAllChannels(sacIdx) & any(strcmp(baseChannels{sacIdx}, allChannels{sIdx}));
    end
end
channels = baseChannels(tfAllChannels); % These are the channels that are common across subjects
baseConditions = allConditions{1};
tfAllConditions = true(1, length(baseConditions));
for sacIdx = 1:length(baseConditions)
    for sIdx = 2:length(subjects)
        tfAllConditions(sacIdx) = tfAllConditions(sacIdx) & any(strcmp(baseConditions{sacIdx}, allConditions{sIdx}));
    end
end
conditions = baseConditions(tfAllConditions); % These are the conditions that are common across subjects

% Create the SELECT variables for channel and condition
for chanIdx = 1:length(channels);
    tfsChannel{chanIdx} = strcmp(dtTimefrequency(:,cChannel), channels{chanIdx});
end
for condIdx = 1:length(conditions);
    tfsCondition{condIdx} = strcmp(dtTimefrequency(:,cCondition), conditions{condIdx});
end

% Display information about subjects
disp('Subjects being merged:');
disp(subjects);

% Loop over conditions
for condIdx = 1:length(conditions)
    condition = conditions{condIdx};
    tfCondition = tfsCondition{condIdx};
    % Loop over channels
    for chanIdx = 1:length(channels)
        channel = channels{chanIdx};
        tfChannel = tfsChannel{chanIdx};
        
        % Display information about where the loop is..
        disp(['        --> ' condition ' / ' channel]);
        
        % Load the data for each subject into memory
        for sIdx = 1:length(subjects);
            tfSubject = tfsSubject{sIdx};
            CurrentRow = tfSubject & tfCondition & tfChannel;
            filenameSource = dtTimefrequency{CurrentRow, cFilename};
            SOURCE(sIdx) = load(filenameSource);
            SPOWER(:,:,sIdx) = SOURCE(sIdx).power;
        end
        
        % In case we have an elektrode, skip it.
        if strcmp(SOURCE(1).channelType, 'Electrode');
            warning('Current channel is an electrode.. skipping.');
            continue;
        end
        
        % AVERAGE OVER SUBJECTS
        power = mean(SPOWER,3);
        % condition and channel are already set..
        subject = newSubject;
        nrOfEpochs = 1;
        warning('Properties such as frequencies, times, channelType are NOT checked.. They are copied from subject source #1.');
        % Write/copy other properties
        channelType =  SOURCE(1).channelType;
        frequencies =  SOURCE(1).frequencies;
        times =  SOURCE(1).times;
        % Write to file and update in database
        filenameDestination = [Database.Path.TimeFrequencyData subject ' ' channel ' ' condition ' ' num2str(min(frequencies)) '-' num2str(max(frequencies)) 'x' num2str(length(times)) '.mat']; %#ok<NOPRT>
        save(filenameDestination, 'power', 'frequencies', 'times',  'channel', 'channelType', 'subject', 'condition', 'nrOfEpochs', '-v7.3');
        UpdateOrAddTimeFrequencyFileInDatabase(filenameDestination, subject, channel, channelType, condition, frequencies, times, nrOfEpochs)
        % Celebrate the succesfulness
        disp(['- Created ' filenameDestination]);
    end
end

end

% 
% 
% 
% 
% 
% 
% 
% % % Select the correct rows
% % tfSubject = strcmp(dtTimefrequency(:,cSubject), sourceSubject);
% % tfCondition = strcmp(dtTimefrequency(:,cCondition), sourceCondition);
% % channels = unique(dtTimefrequency(tfSubject, cChannel));
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% %% BELOW IS THE OLD CODE FOR TFC's...
% 
% 
% 
% function CreateAverageSubject()
% %CREATEAVERAGESUBJECT Get the data per condition and merge the subjects
% %into one average subject
% 
% load('TFData.mat');
% Path = 'D:\ECoG\Data\TimeFrequency Data\';
% %% Get the uniques
% subjects = sort(unique(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Subject'))));
% for sI = 1:length(subjects), tfS{sI} = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Subject')),subjects{sI}); end
% conditions = sort(unique(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition'))));
% for coI = 1:length(conditions), tfCo{coI} = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')),conditions{coI}); end
% channels = sort(unique(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Channel'))));
% for chI = 1:length(channels), tfCh{chI} = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Channel')),channels{chI}); end
% %% Preselect all different types of conditions
% conditions_FaceHouseIdx = false(size(conditions));
% conditions_FaceHouseFHIdx = false(size(conditions));
% conditions_SphereVoluntaryIdx = false(size(conditions));
% for condI = 1:length(conditions), 
%     conditions_FaceHouseIdx(condI) = ...
%         ~isempty(strfind(conditions{condI}, 'FaceHouse'))...
%         & ( isempty(strfind(conditions{condI}, ' Face ')) ...
%         & isempty(strfind(conditions{condI}, ' House ')) );
% % ~eFH & ~(~eF | ~eH)
%     conditions_FaceHouseFHIdx(condI) = ~isempty(strfind(conditions{condI}, 'FaceHouse'));
%     conditions_SphereVoluntaryIdx(condI) = ~isempty(strfind(conditions{condI}, 'Sphere Voluntary'));
% end
%  conditions_SphereIdx = ~conditions_FaceHouseFHIdx & ~conditions_SphereVoluntaryIdx;
% 
% conditions_FaceHouse = false(size(mergedfreqbins,1),1);
% for condI = conditions(conditions_FaceHouseIdx)', conditions_FaceHouse = conditions_FaceHouse | strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')), cell2mat(condI)); end
% conditions_Sphere = false(size(mergedfreqbins,1),1);
% for condI = conditions(conditions_SphereIdx)', conditions_Sphere = conditions_Sphere | strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')), cell2mat(condI)); end
% 
% 
% conditions_AmbiguousIdx = false(size(conditions));
% for condI = 1:length(conditions), conditions_AmbiguousIdx(condI) = ~isempty(strfind(conditions{condI}, 'Ambiguous')); end
% conditions_UnambiguousIdx = false(size(conditions));
% for condI = 1:length(conditions), conditions_UnambiguousIdx(condI) = ~isempty(strfind(conditions{condI}, 'Unambiguous')); end
% conditions_StimulusIdx = false(size(conditions));
% for condI = 1:length(conditions), conditions_StimulusIdx(condI) = ~isempty(strfind(conditions{condI}, 'Stimulus')); end
% conditions_ResponseIdx = false(size(conditions));
% for condI = 1:length(conditions), conditions_ResponseIdx(condI) = ~isempty(strfind(conditions{condI}, 'Response')); end
% conditions_AmbiguousResponseIdx = conditions_AmbiguousIdx & conditions_ResponseIdx;
% conditions_UnambiguousResponseIdx = conditions_UnambiguousIdx & conditions_ResponseIdx;
% conditions_UnambiguousStimulusIdx = conditions_UnambiguousIdx & conditions_StimulusIdx;
% 
% conditions_AmbResp = false(size(mergedfreqbins,1),1);
% for condI = conditions(conditions_AmbiguousResponseIdx)', conditions_AmbResp = conditions_AmbResp | strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')), cell2mat(condI)); end
% conditions_UnambResp = false(size(mergedfreqbins,1),1);
% for condI = conditions(conditions_UnambiguousResponseIdx)', conditions_UnambResp = conditions_UnambResp | strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')), cell2mat(condI)); end
% conditions_UnambStim = false(size(mergedfreqbins,1),1);
% for condI = conditions(conditions_UnambiguousStimulusIdx)', conditions_UnambStim = conditions_UnambStim | strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Condition')), cell2mat(condI)); end
% 
% %% Preselect subjects
% subjects_CW = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Subject')), 'CW');
% subjects_JF = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Subject')), 'JF');
% subjects_DS = strcmp(mergedfreqbins(:,strcmp(mergedfreqbins_header,'Subject')), 'DS');
% 
% %% Start merging the stuff..
% 
% disp('Creating TfcSet ''CWSphere''');
% CWSphere = CreateTfcSetFrom(mergedfreqbins_header,...
%     mergedfreqbins(subjects_CW & conditions_AmbResp & conditions_Sphere,:),...
%     mergedfreqbins(subjects_CW & conditions_UnambResp & conditions_Sphere,:),...
%     mergedfreqbins(subjects_CW & conditions_UnambStim & conditions_Sphere,:),...
%     Path);
% disp('Creating TfcSet ''CWFaceHouse''');
% CWFaceHouse = CreateTfcSetFrom(mergedfreqbins_header,...
%     mergedfreqbins(subjects_CW & conditions_AmbResp & conditions_FaceHouse,:),...
%     mergedfreqbins(subjects_CW & conditions_UnambResp & conditions_FaceHouse,:),...
%     mergedfreqbins(subjects_CW & conditions_UnambStim & conditions_FaceHouse,:),...
%     Path);
% disp('Creating TfcSet ''DSFaceHouse''');
% DSFaceHouse = CreateTfcSetFrom(mergedfreqbins_header,...
%     mergedfreqbins(subjects_DS & conditions_AmbResp & conditions_FaceHouse,:),...
%     mergedfreqbins(subjects_DS & conditions_UnambResp & conditions_FaceHouse,:),...
%     mergedfreqbins(subjects_DS & conditions_UnambStim & conditions_FaceHouse,:),...
%     Path);
% disp('Creating TfcSet ''JFSphere''');
% JFSphere = CreateTfcSetFrom(mergedfreqbins_header,...
%     mergedfreqbins(subjects_JF & conditions_AmbResp & conditions_Sphere,:),...
%     mergedfreqbins(subjects_JF & conditions_UnambResp & conditions_Sphere,:),...
%     mergedfreqbins(subjects_JF & conditions_UnambStim & conditions_Sphere,:),...
%     Path);
% 
% avgCW = Merge(CWSphere, CWFaceHouse);
% avgSphere = Merge(CWSphere, JFSphere);
% avgFaceHouse = Merge(CWFaceHouse, DSFaceHouse);
% avg = Merge(avgCW, DSFaceHouse, JFSphere);
% 
% data = CWSphere;
% data.name = 'CWSphere';
% save('CWSphere.tfcset', 'data', '-v7.3')
% data = CWFaceHouse;
% data.name = 'CWFaceHouse';
% save('CWFaceHouse.tfcset', 'data', '-v7.3')
% data = DSFaceHouse;
% data.name = 'DSFaceHouse';
% save('DSFaceHouse.tfcset', 'data', '-v7.3')
% data = JFSphere;
% data.name = 'JFSphere';
% save('JFSphere.tfcset', 'data', '-v7.3')
% data = avgCW;
% data.name = 'avgCW';
% save('avgCW.tfcset', 'data', '-v7.3')
% data = avgSphere;
% data.name = 'avgSphere';
% save('avgSphere.tfcset', 'data', '-v7.3')
% data = avgFaceHouse;
% data.name = 'avgFaceHouse';
% save('avgFaceHouse.tfcset', 'data', '-v7.3')
% data = avg;
% data.name = 'avg';
% save('avg.tfcset', 'data', '-v7.3')
% 
% end
% 
% function X = CreateTfcSetFrom(RowHeader, AmbRespRows, UnambRespRows, UnambStimRows, Path)
% 
% X.AmbRespData = CreateTfcSingleFrom(RowHeader, AmbRespRows, Path);
% X.UnambRespData = CreateTfcSingleFrom(RowHeader, UnambRespRows, Path);
% X.UnambStimData = CreateTfcSingleFrom(RowHeader, UnambStimRows, Path);
% 
% tfF1 = FrequencyBinsUnequal(X.AmbRespData.FrequencyBins, X.UnambRespData.FrequencyBins);
% tfF2 = FrequencyBinsUnequal(X.AmbRespData.FrequencyBins, X.UnambStimData.FrequencyBins);
% if tfF1 || tfF2
%     error('FrequencyBin mismatch!');
% end
% 
% tfT1 = TimesUnequal(X.AmbRespData.Times, X.UnambRespData.Times);
% tfT2 = TimesUnequal(X.AmbRespData.Times, X.UnambStimData.Times);
% if tfT1 || tfT2
%     error('Times mismatch!');
% end
% 
% end
% 
% %% Create TFC Single from files
% function X = CreateTfcSingleFrom(RowHeader, Rows, Path)
% 
% %for every channel, merge all given conditions
% X.Channels = sort(unique(Rows(:,strcmp(RowHeader,'Channel'))));
% 
% for chI = 1:length(X.Channels)
%     selectedRows = Rows(strcmp(Rows(:,strcmp(RowHeader,'Channel')),X.Channels{chI}),:);
%     for rI = 1:size(selectedRows,1) % Runs per channel
%         DATA(rI) = load([Path selectedRows{rI,strcmp(RowHeader,'Datafile')}]);
%     end
%     
%     if size(selectedRows,1) ~= 1, error('Should be exactly one file for now.. (TODO)'); end
%         
%     % Data voorbereiden
%     isStimulus = isempty(strfind(DATA(rI).condition, 'Response'));
%     [D1 epochSelection] = SelectEpochs(DATA(rI).data, isStimulus);
% %     countS = sum(epochSelection); count = length(epochSelection); countR = count - countS;
% %     disp([subjects{sI} ' ' dset(dsI).cond(condI).condition ' --> Selected ' num2str(countS) ' out of ' num2str(count) ' epoch(s)...  (Removing ' num2str(countR) ')']);
%     D2 = CorrectBaseline(D1, isStimulus);
%     
%     % Frequencybins controleren!
%     if ~any(strcmp(fieldnames(X), 'FrequencyBins'))
%         X.FrequencyBins = DATA(rI).data.frequencybins;
%     else
%         if FrequencyBinsUnequal(X.FrequencyBins, DATA(rI).data.frequencybins)
%             error('FrequencyBin mismatch!');
%         end
%     end
%     
%     % Times controleren
%     if ~any(strcmp(fieldnames(X), 'Times'))
%         X.Times = DATA(rI).data.times;
%     else
%         if TimesUnequal(X.Times, DATA(rI).data.times)
%             error('Times mismatch!');
%         end
%     end
%     
%     % Er is max 1 datafile (zie error 'Should be exactly one file for
%     % now..') dus er hoeft niet te worden geconcat over files.    
%     freqByTimes = mean(D2.powerCorrected,3);
%     X.Power(:,:,chI) = freqByTimes;
%     % Power is a freqbins x times x channels matrix
%         
%     warning('TODO: Implement statistics here!');
% end
% 
% 
% end
% 
% %% Compare
% function tf = FrequencyBinsUnequal(f1, f2)
% tf = any(any([f1.min;f1.max] ~= [f2.min;f2.max]));
% end
% 
% function tf = TimesUnequal(t1, t2)
% tf = any(any(t1 ~= t2));
% end
% 
% function channels = CreateMergedChannelList(ch1, ch2, ch3)
% if exist('ch3', 'var')
%     chAll = cat(1,ch1, ch2, ch3);
% else
%     chAll = cat(1,ch1, ch2);
% end
% chUnique = unique(chAll);
% tf = false(size(chUnique));
% for i = 1:length(chUnique)
%     tf(i) = ...
%         any(strcmp(ch1,chUnique{i})) && any(strcmp(ch2,chUnique{i})) ...
%         && (~exist('ch3','var') || any(strcmp(ch3,chUnique{i})));
% end
% channels = chUnique(tf);
% end
% 
% function idx = CreateChannelIndex(FromChannels, ToChannels)
% idx = nan(size(ToChannels));
% for i = 1:length(ToChannels)
%     idx(i) = find(strcmp(FromChannels, ToChannels{i}));
% end
% end
% 
% %% Merge
% function X = Merge(X1, X2, X3)
% 
% % Make sure the frequency bins are the same
% tfF1 = FrequencyBinsUnequal(X1.AmbRespData.FrequencyBins, X2.AmbRespData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.UnambRespData.FrequencyBins, X2.UnambRespData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.UnambStimData.FrequencyBins, X2.UnambStimData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.AmbRespData.FrequencyBins, X1.UnambRespData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.AmbRespData.FrequencyBins, X1.UnambStimData.FrequencyBins);
% if exist('X3', 'var')
%     tfF2 = FrequencyBinsUnequal(X1.AmbRespData.FrequencyBins, X3.AmbRespData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.UnambRespData.FrequencyBins, X3.UnambRespData.FrequencyBins) ...
%     || FrequencyBinsUnequal(X1.UnambStimData.FrequencyBins, X3.UnambStimData.FrequencyBins);
% else
%     tfF2 = false;
% end
% if tfF1 || tfF2
%     error('FrequencyBin mismatch!');
% else
%     MergedFrequencyBins = X1.AmbRespData.FrequencyBins;
% end
% 
% % Make sure the times are the same
% tfT1 = TimesUnequal(X1.AmbRespData.Times, X2.AmbRespData.Times) ...
%     || TimesUnequal(X1.UnambRespData.Times, X2.UnambRespData.Times) ...
%     || TimesUnequal(X1.UnambStimData.Times, X2.UnambStimData.Times) ...
%     || TimesUnequal(X1.AmbRespData.Times, X1.UnambRespData.Times) ...
%     || TimesUnequal(X1.AmbRespData.Times, X1.UnambStimData.Times);
% if exist('X3', 'var')
%     tfT2 = TimesUnequal(X1.AmbRespData.Times, X3.AmbRespData.Times) ...
%     || TimesUnequal(X1.UnambRespData.Times, X3.UnambRespData.Times) ...
%     || TimesUnequal(X1.UnambStimData.Times, X3.UnambStimData.Times);
% else
%     tfT2 = false;
% end
% if tfT1 || tfT2
%     error('Times mismatch!');
% else
%     MergedTimes = X1.AmbRespData.Times;
% end
% 
% 
% % Make sure we leave out channels that aren't present in EVERY source..
% ch1 = CreateMergedChannelList(X1.AmbRespData.Channels, X1.UnambRespData.Channels, X1.UnambStimData.Channels);
% ch2 = CreateMergedChannelList(X2.AmbRespData.Channels, X2.UnambRespData.Channels, X2.UnambStimData.Channels);
% if exist('X3', 'var')
%     ch3 = CreateMergedChannelList(X3.AmbRespData.Channels, X3.UnambRespData.Channels, X3.UnambStimData.Channels);
%     MergedChannelList = CreateMergedChannelList(ch1, ch2, ch3);
% else
%     MergedChannelList = CreateMergedChannelList(ch1, ch2);
% end
% % Create the indexers to translate data to new channel list
% x1arIdx = CreateChannelIndex(X1.AmbRespData.Channels, MergedChannelList);
% x1urIdx = CreateChannelIndex(X1.UnambRespData.Channels, MergedChannelList);
% x1usIdx = CreateChannelIndex(X1.UnambStimData.Channels, MergedChannelList);
% x2arIdx = CreateChannelIndex(X2.AmbRespData.Channels, MergedChannelList);
% x2urIdx = CreateChannelIndex(X2.UnambRespData.Channels, MergedChannelList);
% x2usIdx = CreateChannelIndex(X2.UnambStimData.Channels, MergedChannelList);
% if exist('X3', 'var')
%     x3arIdx = CreateChannelIndex(X3.AmbRespData.Channels, MergedChannelList);
%     x3urIdx = CreateChannelIndex(X3.UnambRespData.Channels, MergedChannelList);
%     x3usIdx = CreateChannelIndex(X3.UnambStimData.Channels, MergedChannelList);
% end
% 
% % Build the new TfcSet
% X.AmbRespData.Channels = MergedChannelList;
% X.AmbRespData.FrequencyBins = MergedFrequencyBins;
% X.AmbRespData.Times = MergedTimes;
% if exist('X3', 'var')
%     X.AmbRespData.Power = mean( ...
%         cat(4, X1.AmbRespData.Power(:,:,x1arIdx), X2.AmbRespData.Power(:,:,x2arIdx), X3.AmbRespData.Power(:,:,x3arIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% else
%     X.AmbRespData.Power = mean( ...
%         cat(4, X1.AmbRespData.Power(:,:,x1arIdx), X2.AmbRespData.Power(:,:,x2arIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% end
% 
% X.UnambRespData.Channels = MergedChannelList;
% X.UnambRespData.FrequencyBins = MergedFrequencyBins;
% X.UnambRespData.Times = MergedTimes;
% if exist('X3', 'var')
%     X.UnambRespData.Power = mean( ...
%         cat(4, X1.UnambRespData.Power(:,:,x1urIdx), X2.UnambRespData.Power(:,:,x2urIdx), X3.UnambRespData.Power(:,:,x3urIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% else
%     X.UnambRespData.Power = mean( ...
%         cat(4, X1.UnambRespData.Power(:,:,x1urIdx), X2.UnambRespData.Power(:,:,x2urIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% end
% 
% X.UnambStimData.Channels = MergedChannelList;
% X.UnambStimData.FrequencyBins = MergedFrequencyBins;
% X.UnambStimData.Times = MergedTimes;
% if exist('X3', 'var')
%     X.UnambStimData.Power = mean( ...
%         cat(4, X1.UnambStimData.Power(:,:,x1usIdx), X2.UnambStimData.Power(:,:,x2usIdx), X3.UnambStimData.Power(:,:,x3usIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% else
%     X.UnambStimData.Power = mean( ...
%         cat(4, X1.UnambStimData.Power(:,:,x1usIdx), X2.UnambStimData.Power(:,:,x2usIdx)) ...
%         ,4);
%     warning('TODO: Implement statistics here!');
% end
% 
% end

